%FEM_LAME2D  two-dimensional finite element method for the Lame-Problem
%

% Initialisation
% E = 2900; nu = 0.4;
% mu = E/(2*(1+nu)); lambda = E*nu/((1+nu)*(1-2*nu));
int_lambda=1;
int_mu=1;
out_lambda=10;
out_mu=10;
Ntri=size(elem,1);
N=size(node,1);
A = sparse(2*size(node,1),2*size(node,1)); 
b = zeros(2*size(node,1),1);
A = sparse(2*N,2*N); u = sparse(2*N,1); b = sparse(2*N,1);
NTI=size(int_elem,1);
NTO=size(out_elem,1);
% %%%%%%%%%%  Assembly matrix
for j = 1:NTI
  I = 2*int_elem(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0]; 
  A(I,I) = A(I,I) +stima3(node(int_elem(j,:),:),int_lambda,int_mu);   
end
for j = 1:NTO
  I = 2*out_elem(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0]; 
  A(I,I) = A(I,I) +stima3(node(out_elem(j,:),:),out_lambda,out_mu);   
end


% %%%%%%% Volume forces   %%%%%%%
for j = 1:NTI
     I = 2*int_elem(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0];
    fs = int_f(sum(node(int_elem(j,:),:))/3,int_lambda,int_mu)';%三角形重心f(x,in_lamda,in_mu,out_lamda,out_mu)
  b(I) = b(I) +det([1,1,1;node(int_elem(j,:),:)'])*[fs;fs;fs]/6;
    IO = intersect(immersed_inter(:,1),int_elem(j,:));%%%%%%%  检测单元是不是包含immersed point
  if size(IO,2)==2  %%%%%三角单元有两个顶点在immersed interface 上
       I = 2*IO([1,1,2,2])-[1,0,1,0];
       
       
       %%%%%%%%%%   精确积分   %%%%%%%%%%%
       %%%%%%%%%%  梯形公式   %%%%%%%%%%%
     q1 = q(node(IO(1),:))';
     q2 = q(node(IO(2),:))';
      b(I) = b(I)-sqrt((node(IO(1),1) -node(IO(2),1))^2+...
          (node(IO(1),2) -node(IO(2),2))^2)*[q1;q2]/2;
      
%      %%%%%%%%%%  抛物（或 Simpson）积分 %%%%%%%%
%       q1= q(node(IO(1),:))'+1/2*4*q(sum(node(IO(1,:),:))/2)';
%       q2= q(node(IO(2),:))'+1/2*4*q(sum(node(IO(1,:),:))/2)';
%       b(I) = b(I)-sqrt((node(IO(1),1) -node(IO(2),1))^2+...
%           (node(IO(1),2) -node(IO(2),2))^2)*[q1;q2]/6;
%       
    
  end
      
      
      
  
end
for j = 1:NTO
  I = 2*out_elem(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0];
  fs = out_f(sum(node(out_elem(j,:),:))/3,out_lambda,out_mu)';%三角形重心f(x,in_lamda,in_mu,out_lamda,out_mu)
  b(I) = b(I) +det([1,1,1;node(out_elem(j,:),:)'])*[fs;fs;fs]/6;

  
%     IO=intersect(immersed_inter(:,1),out_elem(j,:));%%%%%%%  检测单元是不是包含immersed point
%   if size(IO,2)==2  %%%%%三角单元至少有两个顶点在immersed interface 上
%         q1 = q(node(IO(1),:))';
%         q2 = q(node(IO(2),:))';
%      
%       I = 2*IO([1,1,2,2])-[1,0,1,0];
%    b(I) = b(I)-sqrt((node(IO(1),1) -node(IO(2),1))^2+...
%           (node(IO(1),2) -node(IO(2),2))^2)*[q1;q2]/2;
%   end
      
      %%%%%%%%%    需要处理的 %%%%%%%%%%%%
end 
  


%%%%%%%%%%%%% jump condition  %%%%%%%%%5

% Neumann conditions
% if ~isempty(neumann)
%   n = (coordinates(neumann(:,2),:) -coodrdinates(neumann(:,1),:))*[0,-1;1,0];
%   for j = 1:size(neumann,1);
%     I = 2*neumann(j,[1,1,2,2]) -[1,0,1,0];
%     gm = g(sum(coordinates(neumann(j,:),:))/2, n(j,:)/norm(n(j,:)))';
%     b(I) = b(I) +norm(n(j,:))*[gm;gm]/2;
%   end
% end
u=sparse(2*N,1);
%%%%%%%%%% Dirichlet condition  %%%%%%%%
u(2*Dirichlet-1) = feval('ux_d',node(Dirichlet,:),out_mu) ;
u(2*Dirichlet)   = feval('uy_d',node(Dirichlet,:),out_mu) ;
% u(2*Dirichlet-1) = ux_d(node(Dirichlet,:),out_mu) ;
% u(2*Dirichlet)   = uy_d(node(Dirichlet,:),out_mu) ;

b = b -A * u;


% Calculating the solution
BoundNodes = unique(Dirichlet);
nd=size(BoundNodes,1);%% number of dirichlet nodes 
Bdnodes=zeros(2*nd,1);
Bdnodes(1:2:2*nd)=2*BoundNodes-1;
Bdnodes(2:2:2*nd)=2*BoundNodes;

FreeNodes = setdiff(1:2*N,Bdnodes);
if size(FreeNodes)>0
    B=A(FreeNodes,FreeNodes);
	u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
end
%%%%%%%true  value%%%%%
uu=zeros(2*N,1);
uu(1:2:2*N)=ux_d(node,out_mu);
uu(2:2:2*N)=uy_d(node,out_mu);
er=u-uu;
error1=max(abs(u-uu))
[error,position]=max(abs(u-uu));
pos=floor(position/2);
maxpos=[ceil(pos/Nx),mod(pos,Nx)];
fprintf('the row of the max error point is %d\n', maxpos(1,1));
fprintf('the column of the max error point is %d\n', maxpos(1,2));
% error=u-uu;
%%%%%%%%%%  计算解 %%%

% figure(1)
% 
% trisurf(elem,u(1:2:size(u,1))+node(:,1), ...
%     u(2:2:size(u,1))+node(:,2), ...
%     zeros(size(node,1),1),  'facecolor','interp');
% view(2), axis equal, axis on, title('computational value'),   drawnow;
%%%%%%%%%% 真实解 %%%%%%%
% figure(2)
% trisurf(elem,uu(1:2:size(u,1))+node(:,1), ...
%     uu(2:2:size(u,1))+node(:,2), ...
%     zeros(size(node,1),1),  'facecolor','interp');
% view(2), axis equal, axis on, title('true value')   ,drawnow;
% er=u-uu;
% %%%%%%%%%%  x方向上的误差
% figure(3)
% trisurf(elem,node(:,1), node(:,2), er(1:2:end)',  'facecolor','interp');
% title('error in x direction'),
% %%%%%%%%%%  y方向上的误差
% figure(4)
% trisurf(elem,node(:,1), node(:,2), er(2:2:end)',  'facecolor','interp');
% title('error in y direction')

% 
[AvE,Eps3,AvS,Sigma3] = ...
    avmatrix(node,elem,u,int_lambda,int_mu,...
               out_lambda,out_mu ,int_elem,out_elem);
show(elem,int_elem,out_elem,node,AvS,u,...
      int_lambda,int_mu,out_lambda,out_mu);
  estimate = aposteriori(node,elem,int_elem,out_elem ,...
            AvE,Eps3,AvS,Sigma3,u,int_lambda,int_mu,out_lambda,out_mu)